Ondrej Cernotik
23 February 2015
The evolution of a qubit dispersively coupled to a thermal cavity that is countinuously monitored is given by the stochastic master equation $$ \mathrm{d}\rho = -i\chi[\sigma_z r_\phi, \rho]\mathrm{d}t+\kappa\{(\bar{n}+1)\mathcal{D}[a]+\bar{n}\mathcal{D}[a^\dagger]\}\rho\mathrm{d}t+\sqrt{\frac{\kappa}{2\bar{n}+1}}\mathcal{H}[(\bar{n}+1)a-\bar{n}a^\dagger]\rho\mathrm{d}W, $$ where $\chi$ is the qubit-cavity coupling, $r_\phi = (ae^{i\phi}+a^\dagger e^{-i\phi})/\sqrt{2}$ is a quadrature of the cavity field which couples to an environment with $\bar{n}$ thermal photons at a rate $\kappa$, $\mathcal{D}[O]\rho = O\rho O^\dagger-\frac{1}{2}(O^\dagger O\rho+\rho O^\dagger O)$ is a Lindblad term, and $\mathcal{H}[O]\rho = (O-\langle O\rangle)\rho+\mathrm{H.c.}$ is the measurment term. Since the QuTiP stochastic solver [1,2] treats the measurement operators also as Lindblad operators, i.e., corresponding to decay, it is necessary to rewrite the Lindblad terms in the stochastic master equation using the measurement operator. We have $$ \kappa(\bar{n}+1)\mathcal{D}[a]\rho+\kappa\bar{n}\mathcal{D}[a^\dagger]\rho = \frac{\kappa}{2\bar{n}+1}\mathcal{D}[(bar{n}+1)a-\bar{n}a^\dagger]\rho +\frac{2\kappa\bar{n}(\bar{n}+1)}{2\bar{n}+1}\mathcal{D}[q]\rho $$ with $q = (a+a^\dagger)/\sqrt{2}$ being the amplitude quadrature of the cavity field.
When the cavity evolves on a time scale much faster than the qubit-cavity coupling, $\kappa\gg\chi$, its dynamics can be adiabatically eliminated from the equation of motion. Using Gaussian adiabatic elimination [3], the qubit dynamics is described by the equation $$ \mathrm{d}\rho_q = \frac{2\chi^2}{\kappa}(2\bar{n}+1)\mathcal{D}[\sigma_z]\rho_q\mathrm{d}t+\sqrt{\frac{2\chi^2}{\kappa(2\bar{n}+1)}}\mathcal{H}[-i(2\bar{n}\cos\phi+e^{-i\phi})\sigma_z]\rho_q\mathrm{d}W. $$ We compare this method of adiabatic elimination with results obtained by density operator expansion [4] (this method is, in fact, used only to obtain the measurement term; the decay has been obtained using projection operator method [5] to account for thermal noise), which gives the equation $$ \mathrm{d}\rho_q = \frac{2\chi^2}{\kappa}(2\bar{n}+1)\mathcal{D}[\sigma_z]\rho_q\mathrm{d}t+\sqrt{\frac{2\chi^2}{\kappa}}\mathcal{H}[e^{-i(\phi+\pi/2)}\sigma_z]\rho_q\mathrm{d}W; $$ see Ref. [3] for details.
The elimination methods and the exact solution are quantitatively compared using the trace distance of the two solutions. For a qubit, the trace distance of two density matrices is conveniently expressed using the expecation values of the Pauli operators $$ D(\rho_1, \rho_2) = \frac{1}{2}\sqrt{\sum_{i\in\{x,y,z\}}(\langle\sigma_i^1\rangle-\langle\sigma_i^2\rangle)^2}, $$ with $\langle\sigma_i^j\rangle = \mathrm{tr}(\rho_j\sigma_i)$. We use ensemble averaging to obtain an average trace distance for each elimination method as a function of time; additional time averaging gives a single figure of merit for the performance of the elimination methods.
References
[1] J. R. Johansson, P. D. Nation, and F. Nori, Computer Physics Communications 183, 1760 (2012).
[2] J. R. Johansson, P. D. Nation, and F. Nori, Computer Physics Communications 184, 1234 (2013).
[3] O. Cernotik, D. Vasilyev, and K. Hammerer, arXiv:1503.07484.
[4] A. Doherty and K. Jacobs, Physical Review A 60, 2700 (1999).
[5] C. Gardiner and P. Zoller, Quantum Noise: A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics (Springer, 2004).
In [1]:
%matplotlib inline
In [2]:
from qutip import qeye, destroy, sigmax, sigmay, sigmaz, smesolve, ket2dm, basis, tensor, thermal_dm
from numpy import linspace, cos, sin, pi, sqrt, exp
from matplotlib.pyplot import subplots, plot
In [3]:
def trace_dist(chi, phi, kappa, nbar, Ntraj, NFock, Nstep, Nsub, tmax, rhoq0):
#***** function trace_dist *****
#This function calculates the average trace distance for Gaussian adiabatic
#elimination and density operator expansion as functions of time.
#Input parameters:
# chi: qubit-cavity coupling,
# phi: phase of the field quadrature in the interaction,
# kappa: intrinsic measurement (and decay) rate,
# nbar: thermal occupation,
# Ntraj: number of generated quantum trajectories,
# NFock: Fock number cutoff,
# Nstep: number of steps for evolution,
# Nsub: number of substeps for the stochastic solver,
# tmax: maximum time,
# rhoq0: initial qubit state.
#Returns:
# tlist: list of times,
# D_gae: trace distance for Gaussian adiabatic elimination,
# D_doe: trace distance for density operator expansion.
#***** System operators *****
tlist = linspace(0,tmax,Nstep)
Id2 = qeye(2)
IdN = qeye(NFock) # identity operators on the qubit and the cavit field
a = tensor(Id2,destroy(NFock)) # field annihilation operator
q = (a+a.dag())/sqrt(2.)
p = -1j*(a-a.dag())/sqrt(2.) # field quadratures
sx = tensor(sigmax(),IdN)
sy = tensor(sigmay(),IdN)
sz = tensor(sigmaz(),IdN) # Pauli matrices
#***** Full dynamics *****
H_full = chi*sz*(q*cos(phi)-p*sin(phi)) # Hamiltonian
c_op_full = [sqrt(2*kappa*nbar*(nbar+1)/(2*nbar+1))*q] # jump operators
sc_op_full = [sqrt(kappa/(2*nbar+1))*((nbar+1)*a-nbar*a.dag())] # measurement and jump operators
e_op_full = [sx, sy, sz] # operators whose evolution we are interested in
rho0 = tensor(rhoq0, thermal_dm(NFock, nbar)) # initial state
#***** Gaussian adiabatic elimination *****
H_gae = 0*Id2
c_op_gae = [sqrt(2*chi*chi*(2*nbar+1)/kappa-(2*chi*chi/(kappa*(2*nbar+1)))*(2*nbar*cos(phi)+exp(-1j*phi))*(2*nbar*cos(phi)+exp(1j*phi)))*sigmaz()]
sc_op_gae = [-1j*sqrt(2*chi*chi/(kappa*(2*nbar+1)))*(2*nbar*cos(phi)+exp(-1j*phi))*sigmaz()]
e_op_gae = [sigmax(), sigmay(), sigmaz()]
#***** Density operator expansion *****
H_doe = H_gae
c_op_doe = [sqrt(4*chi*chi*nbar/kappa)*sigmaz()]
sc_op_doe = [sqrt(2*chi*chi/kappa)*exp(-1j*(phi+pi/2.))*sigmaz()]
e_op_doe = e_op_gae
#***** Generating trajectories *****
D_gae = 0
D_doe = 0
for i in xrange(0,Ntraj):
print i
sol_full = smesolve(H_full, rho0, tlist, c_op_full, sc_op_full, e_op_full, nsubsteps=Nsub, method='homodyne',
solver='fast-milstein', store_measurement=True)
sol_gae = smesolve(H_gae, rhoq0, tlist, c_op_gae, sc_op_gae, e_op_gae, nsubsteps=Nsub, method='homodyne',
solver='fast-milstein', noise=sol_full.noise)
sol_doe = smesolve(H_doe, rhoq0, tlist, c_op_doe, sc_op_doe, e_op_doe, nsubsteps=Nsub, method='homodyne',
solver='fast-milstein', noise=sol_full.noise)
s1x = sol_full.expect[0]
s1y = sol_full.expect[1]
s1z = sol_full.expect[2]
s2x = sol_gae.expect[0]
s2y = sol_gae.expect[1]
s2z = sol_gae.expect[2]
s3x = sol_doe.expect[0]
s3y = sol_doe.expect[1]
s3z = sol_doe.expect[2]
D_gae = D_gae+sqrt((s1x-s2x)**2+(s1y-s2y)**2+(s1z-s2z)**2)/2.
D_doe = D_doe+sqrt((s1x-s3x)**2+(s1y-s3y)**2+(s1z-s3z)**2)/2.
D_gae = D_gae/Ntraj
D_doe = D_doe/Ntraj
return tlist, D_gae, D_doe
In [4]:
chi = 0.1
phi = pi/2.
kappa = 1.
nbar = 0.5
Ntraj = 20
NFock = 8
Nsteps = 1000
Nsub = 250
tmax = 50
t, D_gae, D_doe = trace_dist(chi, phi, kappa, nbar, Ntraj, NFock, Nsteps, Nsub, tmax, ket2dm((basis(2,0)+basis(2,1)).unit()))
In [5]:
fig, axes = subplots(figsize=(12,8))
axes.plot(t, D_gae, label='Gaussian ad. elimination')
axes.plot(t, D_doe, label='Density op. expansion')
axes.set_xlabel('Time')
axes.set_ylabel('Trace distance')
axes.legend(loc=0)
Out[5]:
In [6]:
D_gae_time = D_gae.sum()/Nsteps
print D_gae_time
In [7]:
D_doe_time = D_doe.sum()/Nsteps
print D_doe_time